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Abstract. This paper presents the algorithm for determining the Lemaitre-Tolman 
model that best fits given datasets for maximum stellar ages, and SNIa luminosities, 
both as functions of redshift. It then applies it to current cosmological data. Special 
attention must be given to the handling of the origin, and the region of the maximum 
diameter distances. As with a previous combination of datasets (galaxy number counts 
and luminosity distances versus redshift), there are relationships that must hold at the 
region of the maximum diameter distance, which are unlikely to be obeyed exactly by 
real data. We show how to make corrections that enable a self-consistent solution to 
be found. We address the questions of the best way to approximate discrete data with 
smooth functions, and how to estimate the uncertainties of the output — the 3 free 
functions that determine a specific Lemaitre-Tolmanmetric. While current data does 
not permit any confidence in our results, we show that the method works well, and 
reasonable Lemaitre-Tolman models do fit with or without a cosmological constant. 



1. Introduction 

Over the past century, the investigation of the structure and evolution of the universe has 
become a major scientific thrust. Einstein's equations made the geometry of the cosmos 
inseparable from its matter content. Thus there are many intriguing investigations, 
geometric as well as physical, ranging from whether the spatial section of universe is 
open or closed, and what the global topology is, down to lensing by individual mass 
concentrations, and detecting the signatures of galactic black holes. Equally important 
is the geometry on a range of intermediate scales, and its effect on the paths of the light 
rays that constitute the bulk of our cosmological observations. 

The assumption of a homogeneous Friedmann-Lemaitre- Robertson- Walker (FLRW) 
universe made the Einstein equations tractable, and later, complemented by 
perturbation theory, opened the way to our present day understanding of an observable 
universe that is 10^^ times larger than a human frame, and 10^ times older than a human 
lifespan. Currently, in most cases the cosmological observations are anlysyed within the 
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framework of homogenous models, for example supernova observations [8] galaxy redshift 
surveys [86] in particular baryon acoutiong oscilations [12], and cosmic microvawve 
backrogund radition [49]. The cosmological principle remains a core ingredient of the 
standard concordance model, but as our measurements become ever more complete 
and detailed, deviations from homogeneity will have to be taken seriously. Certainly, 
inhomogeneity exists on a wide range of scales, and a knowledge of the associated 
spacetime curvature is indispensable to the detailed analysis of observational data. 

Whether one views it as detecting what level of inhomogeneity exists, or verifying 
how well homogeneity is satisfied, the investigation is essential. The most important 
reason is that the scientific method requires it — if you can test an assumption you 
should — not just in one way, but in many different ways. 

The history of science can provide many examples of why different tests are 
important. The most well known example is probably the Ptolemaic vs Copernican 
model of the Universe. At the time of Nicolaus Copernicus, the Ptolemaic model with 
its large number of epicycles was able to fit cosmological/ astronomical data much better 
than the Copernican model. Thus, we should bear in mind that consistency with current 
observations is not a final proof oi the underlying assumptions of the model, and as more 
data become available we should test our assumptions using different methods. 

One way of testing the assumption of homogeneity is to look at the consistency 
relations between the observables. A consistency check is based on verifying whether 
the relation between observations is as given by the homogeneous cosmological model. 
Previously, Ribeiro and Stocgcr considered the consistency between the expected 
number counts in a homogeneous model, a commonly used galaxy luminosity function, 
and the observed galaxy number counts [74]. Later they showed that such an analysis 
strongly depends on the distance definition used [3]. In [25] Clarkson, Bassett & Lu 
showed that in Priedmann models that H{z) and D{z) are not independent but must 
obey a special relation. Thus, by checking if observations obey this relation, one can 
test the assumption of homogeneity. A limited version of such a test was proposed in 
[78], where just H{z) alone can be used to test the consistency between observations 
and a spatially flat homogeneous model. Another interesting aspect is to check the 
consistency between the age of the universe inferred from the cosmological models with 
'local' (age of meteorites or stars in our own Galaxy) measurements. 

In this paper we take another approach. We study the observations within an 
inhomogeneous framework. The motivation is as follows: if the the Universe is indeed 
homogeneous, then the data should favour a homogeneous model over inhomogeneous 
ones. 

As argued in [16], spherically symmetric models are the logical flrst step away from 
homogeneity. Radial inhomogeneity is intrinsically hard to detect and separate from 
other past null cone effects. Our observations reveal a tilted slice through spacetime, and 
what we see is affected by the details of the cosmic evolution due to the bulk equations 
of state, the evolution of light sources, and the lumpiness of the matter distribution. 

The simplest spherically symmetric, non-stationary, non-homogeneous solution of 
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the Einstein equations is the Lemaitre-Tolman (LT) model. Currently there is a large 
number of papers that use the LT model to fit cosmological data without dark energy 
— a non-compete list includes [22, 48, 7, 24, 6, 36, 10, 20, 21, 59, 13, 58, 19, 39, 40, 
28, 29, 35, 11, 95, 17, 37, 4, 88, 73, 26, 64, 54, 94, 92, 93, 67, 32]; for a review see 
[15, 60]. However, these methods are based on what is often called the direct approach, 
which can be summarized as follows: 'assume a model, fit the model to the data, 
obtain the best fit parameters of the model and compare the goodness of fit with the 
standard cosmological model'. The caveat of this kind of approach is that any assumed 
parametrization is not general, and may omit a large part of the set of all possible 
models, hence a better approach is the so called inverse approach. Within the inverse 
approach one uses the data to construct the model. In theory, this approach is much 
less restricted by choice of model and parameterization than the direct approach. In 
practice however one still needs to impose some additional assumptions to make the 
calculations tractable. The problem is that we need to be able to solve the Einstein 
equations. Thus, one either employs a series expansion around the observer [53, 68, 69], 
or assumes spherical symmetry to calculate the evolution of the Universe within the 
whole observable domain. Here one can either use 'observer coordinates' and the fluid- 
ray tetrad approach [34, 84, 57, 9, 45], or use standard coordinates in the spherically 
symmetric dust model, i.e. the LT model, and determine the null-cone data relations 
[66, 43, 56, 63, 91, 23, 90, 75, 76, 77]. 

In this paper we will use the inverse method with the Lemaitre-Tolman model to 
test the assumption of homogeneity. The motivation and plan is as follows: 'start with 
an inhomogeneous model more general than FLRW, apply the inverse method to the 
actual observational data to specify the model, then if the model obtained in this way 
is homogeneous (or very close to homogeneous) we find in favour of homogeneity, if on 
the other hand the model is significantly inhomogeneous, we disfavor the homogeneity'. 
It should be emphasised that we can only favour or disfavor, not actually prove or 
rule out. As shown in Ref. [66] a Lemaitre-Tolman model could be found to fit any 
reasonable observational data for the luminosity (or diameter) distance versus redshift, 
and the galaxy number-counts versus redshift; while the difficulty of separating cosmic 
evolution, source evolution, and inhomogeneity was pointed out. This concept was 
expanded, with an important correction, and explicit numerical implementations were 
developed in [43, 56, 63]. A possible method for testing source evolution theories in the 
presence of inhomogeneity was suggested in [41]. 

It is important to begin making the connection with real data, as we do here. We 
note however, that the quantity and quality of current data is not sufficient to draw any 
clear conclusions from the LT models obtained. The purpose is to explore the various 
options and uncover potential difficulties, so that a proper analysis of future data will 
be possible. 

It is also important to investigate the question of homogeneity using many different 
methods and datasets. Consequently, this paper applies the metric of the cosmos 
approach to a different set of observations. It uses the luminosity distances of supernovae 
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and the maximum ages of galaxies, both as functions of redshift. This is in contrast 
to using the redshift variations of luminosity or diameter distances and source number 
counts, as considered previously [42, 56, 63]. 



2. The Lemaitre-Tolman solution 

The LT (Lemaitre-Tolman) model [55, 87] is a spherically symmetric solution of 
Einstein's equations where the gravitational source is dust. In comoving and 
synchronous coordinates, the metric is 

ds^ = -dt'' + y^^dr^ + R\t. r){de'' + sin^ ed<p') , (1) 

where W{r) = 1 + /(r) is an arbitrary function, R' = dR/dr. The spatial sections are 
flat only if W{r) = 1, so FT is a geometric factor. The areal radius R{t,r) obeys 

R'^W'-l + ^ + ^R\ (2) 
R o 

where R = dR/dt and A is the cosmological constant. From (2), E = [W'^ — l)/2 is the 
total energy per unit mass of the dust particles. Equation (2) is a first integral of the 
Einstein equations, and M(r) is another arbitrary function of integration that gives the 
gravitational mass within each comoving shell of coordinate radius r. The mass density 
in energy units is: 

2M' 

where k — SttG/c^. Equation (2) can be solved by simple integration giving a relation 
for the age of the Universe, r 

(4) 



r{r)^t-tB{r)= f 



- 1 + 2M/R + |Ai?2 

where tB(r') appears as an integration function, and is an arbitrary function of r giving 
the time at which i? = on each worldline. This means that the Big Bang is not 
simultaneous as in the FLRW models, but occurs at different times at different distances 
from the origin. Since all the formulae given so far are covariant under coordinate 
transformations of the form r = g{r), one of the functions W{r), M{r) and ^^(r) can 
be fixed at will by the choice of gf. J Therefore, once this choice is made, a given L-T 
model is fully determined by two of these arbitrary functions. In other words, we need 
two observed relationships to pin down an LT model. For more information on the LT 
model, see [51, 71, 44]. 

I But note that no one choice can cover all possible LT models; any given M{r) either does or does 
not have a vacuum region, it cither does or does not have a spatial maximum; any /(r) either does or 
does not change sign; any tsir) either does or does not have a constant region. 
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2.1. Past Null Cone and Observations 



— = 1 



(6) 



Since light propagates on null geodesies, let us define the past null cone (PNC) by the 
observation event {t, r) = {to, 0) for a central observer at the present day.§ The equation 
for incoming radial null geodesies is 

dr~ W ■ 

Let us choose the radial coordinate so that, on the PNC, 

R 
W 

where by a hat or a subscript wedge (e.g. Q or Qa) we denote any quantity evaluated on 
the PNC. (Note: this choice of r is possible only on a single light cone. In the following, 
we always refer to the PNC of {to, 0).) This gauge choice simphfies the PNC equation: 

?(r) = to-r . (7) 

The equation for the redshift reads [18, 66] 

1 Mr) ^ ^ . . 

1 + z{r) dr W ' ^ ^ 

We aim to show how to define an LT model from luminosity and age data, so we 
need to express the null cone relations using observables like luminosity (or diameter) 
distance and the age of the Universe. In the LT model with observer at the origin, the 
angular diameter distance equals the areal radius R at the point of emission. Using the 
reciprocity theorem [38, 70, 33] this can be related to the luminosity distance: 

Dl 



R = Da = , .2 • 

(1 + z)'^ 

Taking the total derivative of the areal radius R gives [66, 44] 

dr 





dR 


and 


dr 


d^^ _ 


\^ w 


dr^ 



W-R 



R" 



i?' 1 - TT^ 




M'R 

wm 



(9) 



(10) 



(11) 



Note that there is a sign error in eq. (3.19) of [44]. However, what we obtain from 
observations are R{z) and t{z), so we thus convert the above to derivatives with respect 
to redshiftz, using 



dr dz (f if 



d^R d^R 1 



dr^ dz'^ ip^ 



dR 1 d^p 
dz (f^ dz 



(12) 
(13) 



§ As stated elsewhere, there axe good reasons for working initially with a central observer: much 
present-day data is averaged over the whole sky just to get a big enough dataset, we do live at the 
centre of our own past null cone, and radial inhomogencity is much harder to detect than anisotropy. 
There is no implicit claim we live at the centre of a gigantic inhomogeneity. 
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and the redshift derivative of 93 is then 

dz + dl I ' ^^^^ 

which can be expressed as 



dz2 



do? 

dz ^ / . 



+ ^ ifB (f M Ar\ ^ ^ 



1 + z) r^\\r2 3 ; " (1 + ^) 

dz ' ^ ' 



(16) 



where 



RR^ Jo R"^ 

The second set of observations, the age measurements are given in the form of t{z): 
r(z) = t(z) - tB(z) ^to- r(z) - tB(r(z)) . (18) 
The derivatives of r with respect to redshift are 

'^^if.^JLA-w'^B. (19) 

dz ^ dz dz dz 

Solving (10) and (2) for W we get 

dr) 2(f) 
which in terms of redshift is 

Differentiating the above with respect to redshift, and using (19) to ehminate dW/dz, 
we get: 



dM -'f.+i'£ + ^WBR--^ 
dz A ^ 



(22) 



R 

The arbitrary LT functions that are determined by the given r{z) and Dl{z) 
(or equivalently by R{z)), can be calculated in the following way. We first solve 
(22) for M{z) [where R is given by (2), R = -M/R? + ^A/3, and ^ from (12) is 
{dR/dz)/{W - i?)] to find M{z) for the next step, i.e. M^+i. Then Wi+i follows from 
(4). Next r{z) is evaluated by solving (14), and then tB^z) follows from (18). It is 
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convenient to let tsiO) — so that to — ''"(O)- While the above are the basis for a 
numerical algorithm over most of the z range, two loci require special attention, the 
origin, and the maximum in R{z) [56]. 

2.2. Origin 

At the origin, R{t, 0) = V t, and similarly R{t, 0) = = R{t, 0) etc. Regularity [46, 65] 
requires M ~ i?"^ and / ~ i?^ as r — > (where / = W'^ — 1 = 2E). It is then evident 
that — 7> 1 and from (6) and (10) that i?' — )> 1 and di?/dr -> 1 here, meaning i? ~ r. 

Equation (2) shows that i? ~ r and R' ~ r", so it follows from (8) that d2;/dr ~ R' . 
Thus, with z = here, several quantities in the above equations go to zero and we have 
0/0, which makes a direct numerical solution difficult at this locus. In order to find the 
solution in the vicinity of the origin, we use a Taylor series approach, and we write 



f^J2^,z\ ^ = T^J2^iz\ M = J]M,zS / = (23) 

i=l i=l i=0 i=3 i=2 

and of course = dr/dz etc. Since there is a well-defined mapping between r &l z near 
the observer, we write 



R^zU , and R = zU . 



(24) 



so that 



Jo 



dU 



2M I / I AU"^ 
z2 ^ 3 



Inserting (23) and re-expanding the integrand as a Taylor series gives| 



IP 



U 2j 



Ma h 
U 2 



2U^ 



z' + 



dU 



(25) 



-77- > (26) 



where the value oi U at z 



is 



2M3 ^ „ AC/' 

+ /2 + 



U 3 
Taking the z derivative of (26) gives 

dr 
dl 



u 2 







u 2 



M4 

u 2) 



M5 U 
u 2 



r/2 



IP 



fs 

U 2 



U 2 



2U^ 



2U^ 



(27) 




z + 



II Of course, (26) does not constitute a Taylor series in z for r, since the upper limit of each integral 
itself depends on z. Rather, (31) and (32) constitute the first two coefficients of the Taylor series for r. 
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At 2; = 0, these become 







R 


Riz 


u 











z 


z 











dU 


\dR U 




z dz z^ 















f^' dU 


To 







































Ri , 

{Ri + 2R2Z) {Riz + R2Z'') 



— R2 , 



fs\dU 

u 2) ill 



From (12) we see that 
h^Ri- 



(29) 

(30) 

(31) 
(32) 

(33) 



We next substitute these series into (15), re-expand them as Taylor series, and solve for 
each power of z in turn. The /j then follow from (21). We find 



^2 



^1 ^ 



1 - 



2M, 



ARj 



2R2 - M3 



i?2M3 ^ 3i?3 ^ g ^ M4 5M3 

r4 = -^3 \ — h -R4 + H 3— (,34) 



2i?i 



8 



2M4 



i?i 3 ~ ' yRi ) Ri 

Substituting for /2 from (35) into (31) with (27), the expression 

^1 dU 



To = 



2M3 _ 2M3 I AC^ _ , 
(7 ^ 3 3 



(35) 



(36) 



allows one to iteratively solve for M3 using known tq and Ri values. Similarly, 
substituting for /3 from (35) into (32) allows M4 to be found via another iterative 
process. Knowing the Mj, one obtains the fi and fj from (35) & (34). Obviously, one 
gets tB from (18) 

(37) 



^S,0 — ^0 ~ To , 



r,- + r,- 



Thus we can evaluate all coefficients at a small distance away from the origin, so that 
the numerical integration can proceed. 



2.3. Maximum 

A common feature of many big-bang cosmologies is that the diameter distance R{z) 
has a maximum Rm that occurs where our past null cone crosses the apparent horizon 
[52, 43, 56, 63, 44]. Subscript m will indicate evaluation at this maximum. By definition, 
this locus has 



dR 
d7 



0, 



(38) 
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and it follows from (10) & (2) that 



Rrr^^Wrr 



2M„ 



while (21), (11) and (22) show 



(PR 



and 



dr^ 



dM 



dz 



M' ' 
WR. 

dr 



d^R 




(p dM' 




dz^ 


m 


WR . 


m 



WB 



M AR\ W 



i?2 



l + z) 













R, 



(39) 



(40) 



.(41) 



Now the If DE (16) has dR/dz in the denominator, but (40) and (41) ensure that the 
corresponding numerator. 



d^R 
dz^ 



(aR - 



1 dr 



M AR 



w 



W dz 



i?2 



(1 + ^) 



0, 



(42) 



is zero. Similarly, (21) contains 0/0, and it is this expression that is understood wherever 
W appears in (22) and (16). So although dc/p/dz and dMjdz are not actually divergent, 
their DEs are not suitable for numerical integration here. Therefore we seek a Taylor 
series solution in the neighbourhood of Rm-, in powers of lS.z — z — Zm, 



i=2 



i=l 



M = M^ + J2MiAz\ W = W^ + J2Wi^^\ f = f^ + J2fAz' 

i=l i=l i=l 

Again expressing the r integrand as a Taylor series gives 



(43) 



Jo \ V2 RjRl 



p2 

1 L-^ 



h , MA 3 //, , M, 



2+ R 



2+ 



R-^ 



dR 

Rm. 



where 



Rm 



2M^ A^ 



R 



(44) 



(45) 



and thus its z derivative is 



The Metric of the Cosmos from Luminosity and Age Data 

where dR/<\z\m = 0, so that the values at the maximum are 

"■^ dR 



-I 

Jo 



Rr, 

Rm 



2 Rj Rl 



The scries solution of the (p DE, multiplied by WdR/dz, produces^ 

—2WmRmR2 



ri 



Ml 
fi 



+ -n—] TTT-r \ RmR2 , 



[1 + Zm) Ml ; 2Mi 



'^S = ^ ^ ITT - 777:r^ T - 



3Mi 



Ml 2(1 + Zm) 



+ 



mo. 1 2M3 2M| 



+ 



+ 



3Mi(l + z„) 2(l + z^)2 Ml Ml 



Ml 



9Mi 



,(1 + -Zm) Ml 

and feeding these into (21) times di?/d2; gives 



1 _!_/?! RmR'2 



?,Wm mil 3Mi ' 



/l 

/2 



1 '.^2 ,,^2 12 



;i + ..„) + ^^''-"''Mi 

(1 + itK' - 2ARl) ml 
2Mf 



2Mi 



2 {kRl-l)M, 
^ Ml 



R2 



R-m 



^ {ARl-l)3Rs ^ {AARl-A)R, ^ 



2Mi 



-Rm(l + Zm) R. 

Now (39) gives M^, and this plus (47) with (45), 



Mi(l + Z^) {l + Zmf 



2Rr, 



Rr, 



dR 



2M„ If, Mi 



R 



allow fm to be found recursively. Similarly, putting /i from (53) into (48) 

I 1 , . R2 \ o Ml Ml 



-/ 



dR 

J- i^rn 



% Note that (49) also satisfies (40) evaluated at Zn 
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permits a recursive solution for Mi, so that /i and fi follow from (53) and (49). 
This applies to all higher orders, but the coefficient expressions rapidly become large. 
However, for ts, (18) leads to 

tB,m ^to- {Vm + Tm) , tB,i = "(r^ + Tj) , (56) 

and is not determined above. It must be fixed by making the maximum-series curve 
for tB{z) join up with the numerical curved for tsiz) a small z-distance before the 
maximum."^ 

These properties of the maximum in R are not merely difficulties to be overcome, 
they provide independent information with which the results of numerical integrations 
may be cross-checked [43, 56, 63]. (In fact (39) applies also in the case of non-zero 
pressure [5].) If one can measure Rm, then, for given A, (39) immediately gives the total 
gravitational mass within this radius. This holds regardless of the magnitude of any 
intervening radial inhomogeneity, and it holds only at Rm- As we will see, this M must 
agree with the M (and W) obtained from numerical integration. 

3. Data Functions 

The algorithm outlined above requires us to provide r, R, d'^R/dz^, etc at any 
given z. Real observations merely give us a set of discrete values. Hence there is a 
problem of how to calculate higher derivatives based on discrete values. In the case of 
ideal measurements one could use numerical algorithms to estimate the derivatives, for 
example the three-point or five-point derivation rules [72]. In practice however, the data 
are too noisy (supernova data at z 0.4 — 1) or there are too few measurements (age 
data, and supernova data for z > 1). 

One method for extracting a smooth curve from the data would be to apply a 
procedure like the moving average. Such a procedure was outlined in Ref. [63]. Another 
would be to apply some kind of a smoothing, like Gaussian smoothing - for details see 
[81]. The question is how to smooth the data without distorting it, or in other words 
how to distinguish between noise and real variability of the data. In the literature this 
problem is mostly investigated in the context of a reconstruction the evolution of dark 
energy where the data also needs to be differentiated twice to constrain the dark energy 
equation of state [1, 30, 2, 89, 31, 80, 81, 27, 47]. 

To highlight the problem, let us for a moment focus on one procedure of smoothing 
the data. Let us consider the supernova data, Dl{z), in a given redshift bin, and apply 
the least squares method to fit a polynomial to the points in the bin. Now, the bin 
is moving, i.e. for each supernova we consider n — 1 supernovae around it and we fit 
a polynomial to these n points. Then we move to the next supernova and repeat the 

With the distance- number counts-redshift method of ear her papers, Wm and Mi could only be 
determined by matching the maximum series to the numerical integration at some point just before the 
maximum. The difference is probably because a differential equation for dM/dz in that method has 
been replaced by an integral equation for r here. However, each M coefficient now requires an iterated 
integration, which was not necessary previously. 
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procedure*. On the one hand, the smaller n is, the smaller the redshift interval, and 
therefore the more the flexibility of a fit — one can extract more detailed information 
about the curve. On the other hand, if n is too small then the fitting procedure picks up 
the 'noise' of the data. This is illustrated in Fig. 1. If the number of supernovae within 
an interval is n = 10, then the noise becomes significant. Even with n = 200 we can 
still see some artifacts like the minimum around z — 1.3, which is rather unphysical. 
This shows that we still need a much larger large number of observations to generate a 
reliable average Dl{z) curve. ji At this point however, when we only have a few hundreds 
of supernova measurements, and a mere 30-odd age estimations. The present data is 
simply not good enough to be used in the algorithm outlined above. The situation will 
obviously improve over the coming years. 

Still, we would like to use the real data to see what kind of results can be obtained 
when applying our algorithm. We therefore decided to apply a crude approach and to 
fit a chosen curve to the whole data. Of course, this is something that one would prefer 
to avoid. The idea of the inverse problem is to use pure data, without any assumptions 
of its form and use it to specify the geometry of the spacetime. Nevertheless, let us 
proceed in this way, bearing in mind that once the data is of better quality and large 
quantity, it will be possible to apply the method in its true spirit. 

Inevitably, the choice of a fitting function must affect the form of the output. In 
order to minimise this bias we consider several different functions, so that we can assess 
the impact of such choices on the results of our analyses. 

• Age data 

The maximum stellar age data we used comes from [62, 82, 83], and we considered 
3 fitting functions: 

(i) Type 1 

t{z) ^to + tiz + t-2Z^ , (57) 

where to = (11.3355 ± 0.7804) x lO'^ y, t^ = (-9.80635 ± 1.559) x 10*^ y, and 
t2 = (2.3029 ± 0.711) X lO'^ y (see Fig. 2). 

(ii) Type 2 

t{z) = ioe-*^^ , (58) 

where t^ = (12.3137 ± 0.9151) x 10^ y and ii = 1.21623 ± 0.1046. 
(ni) Type 3 

where a = (14.227 ± 1.318) x 10^ y, and h = 2.05036 ± 0.1873. The 
above form is suggested by the Einstein-de Sitter time behaviour, where 
T = (2/3)//o(l + ^)-'/'. 

* The size of the interval will thus vary with redshift. Alternatively one can fix the redshift interval, 
but then the number of points in each interval will no longer be constant. 

ft To investigate angular variation, or equivalently check isotropy, we would need far more data in each 
directional bin. 
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□ 



Q 




Figure 1. Supernova data (upper left panel) and the different fits (other panels). 
For each graph, the value of n denotes the bin size: within each redshift interval that 
contains n supernovae, a 3rd order polynomial is fitted to the data. 



• Angular diameter distance data 



The luminosity data, usually given in the form of distance moduli, /x, where 

Da 

l^^^^og + 25 , 

(1 + z)^ 

come from the supernova observations, here the Union data set [50]. For our 
purpose, it is more convenient to express the functions in the form of the angular 
diameter distance Da{z). We consider 5 different forms of fitting function: 

(i) Type 1 

Da{z) = zDic-^^' , (60) 

where Di = 4.13529 ± 0.1195 Gpc and D2 = 0.909403 ± 0.03956 (see Fig. 2). 

The Hubble constant for this fit is 

-1 

= 72.496 ± 2.095 km s"^ Mpc" 

0/ 



Ho 

(ii) Type 2 



dD, 



dz 



Da{z) = Diz + D^z^ + Dsz^ , (61) 

where Di = 4.04282 ± 0.1534 Gpc, D2 = -3.2535 ± 0.3655 Gpc, and 
D2 = 0.863957 ± 0.2028 Gpc. The Hubble constant for this fit is Hq = 
74.154 ± 2.814 km s"^ Mpc'^ 
(iii) Type 3 

Da = atan (zB) , (62) 

1 + z 
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Figure 2. Left: The maximum stellar age data, and the best-fit polynomial (quadratic) 
of eq. (57), one of 3 different fitting functions tried. Allowing for star formation time, 
the dashed line is the best fit + 1.21 Gy. Right: Angular diameter distance from 
supernova data, the solid line represents the best-fit relation (60), the first of 5 functions 
tried. 



where A = 4.18101 ± 0.293 Gpc, B = -1.0185 ± 0.0913. The Hubble constant 
for this fit is Hq = 70.401 ± 8.822 km s~^ Mpc"^ 

(iv) lype 4 

Da^Ax atan (zB) , (63) 

where A = 1.34115 ± 0.04135 Gpc and B = 2.95013 ± 0.2161. The Hubble 
constant for this fit is Hq = 75.771 ± 6.022 km s~^ Mpc~^ Note that the 
function atan(a;) has a maximum at infinity, so no correction for systematics 
(see below) can be apphed for these 2 choices. As long as the data do not 
approach Rm, there should be no problem. 

(v) Type 5 

Da = Dkcdm , (64) 
where D\cdm is the angular diameter distance as in the standard cosmological 
model - the isotropic and homogeneous ACDM model. We use this case to 
check how this choice (which is the standard choice) influences the analysis. 



4. Correction for systematics 

It was shown in [43, 63] that the functions from the maximum series solution will not 
match properly with the nearby numerical solutions, if the data contains any systematic 
error. In other words, the data as measured does not correspond to a fully self-consistent 
metric solution, in the sense that the various relations that must hold at will not 
be exactly true at Zm, or they will not all hold at quite the same z. This will be the 
case with all real data. It was further demonstrated how these area-distance-maximum 
relations can be used to estimate a "correction" to the data so that a consistent solution is 
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obtained. The same consistency requirement arises for the present datasets and solution 
algorithm; they are the second of (39) and (42). 

If we solve the equations without correction, then these relations will not be exactly 
satisfied, and we will encounter an undesirable numerical behaviour near the maximum 
— see Fig. 3. The point is that at Rm several relations must hold, and imperfect data 
will not satisfy them. The observational data is always accompanied by systematics, 
therefore a correction is needed if a vahd solution is to be obtained. Note that this check 
detects the presence of a cumulative systematic error at Rm, and enables one to find 
a correction that makes the data self-consistent. However, the Rm conditions do not 
show what correction is needed away from the maximum, so an appropriate smoothness 
assumption must be made. 

Since the age data is probably the most sensitive to systematics. For this paper, we 
boldly assume that the age of the oldest stars at each radius tracks the local age of the 
Universe. In fact, the oldest observed stars, most probably, are not even of Population 
III, which are believed to be very massive and rapidly evolving. Therefore we focus on 
"correcting" the age data, and we propose the simplest correction, i.e. adding a constant 
to the age data, which corresponds to the fact that the oldest stars observed in galaxies 
formed some time after the Big Bang. 

To demonstrate how the correction procedure is working, let us consider a model 
with t{z) and Da{z) data of type 1, i.e. as given by (57) and (60) respectively, and a 
cosmological constant, A, corresponding to Q\ = 0.7. 

Fig. 3 presents the results obtained without our correction procedure. The Taylor 
expansion solution in the vicinity of the maximum is plotted using a dashed line, the 
numerical integration with a solid line. As seen, the two results do not match smoothly. 
However, correcting the age data by adding a constant, i.e. r — > r + St, results in a 
smooth fit. The functions M{z) and W{z) after the correction are presented in Fig. 4. 
The correction (St) turns out to depend on the value of A. For observational data in form 
of (57) and (60) with = 0.7 (where Qa = {c^A)/i3H^)), we find 6r = 1.21 x 10^ y; 
for Qa = 0.4 we find 6r = 0.94 x lO'' y, and for Qa = 0, 6t = 0.65 x 10^ y. Given 
the scatter in the data, and likely nature of the observed oldest stars, these values are 
surprisingly reasonable. 

5. Algorithm 

The above was turned into a numerical procedure, written in FORTRAN. The following 
gives the key elements of the algorithm used. 

(i) In the neighbourhood of the origin, fit low order polynomials to Da{z) and r{z), and 
solve for the Taylor coefficients in (34)-(37), using a bisection method, as described 
at the end of section 2.2. These provide initial conditions at a small z value, for 
the numerical integration of the DEs, which can't be evaluated ai z = 0. 

(ii) Choose functional forms for Dl[z) and t[z), and do least-squares fits to the 
observational data. 
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Figure 3. Results obtained witti uncorrected data, for a model with = 0.7. The 
dashed line is the series solution near the maximum R^, and the solid line is the 
numerical solution. It is evident that the series and numerical solutions do not match 
smoothly. Left: M{z). Right: W{z). 




Figure 4. Results after correcting the age data, showing the numerical integration 

(solid line) together with the series solution near Rm (dashed line) for the n\ = 0.7 
case, where the correction to the age data is of form to ^0+^''') where At = 1.21 x 10^ 
y. Left: M{z). Right: W{z). 



(iii) At each step (each zi) of the numerical integration, solve (22) to find the value of 

M for the next step, Mj+i. 

(iv) Use Mi+i and (4) to find Wi+i. 

(v) Solve (14) to find rj+i. 

(vi) Obtain tsi from (18). 

(vii) In practice steps iii-vi can be combined into one integration routine, calculating all 
function values at each z step. We did this by means of a 4-th order Runge-Kutta 
method. 

(viii) In the neighbourhood of the maximum, i?^, do a low order polynomial fit to Dl{z) 
and t{z), and then calculate the coefficients of the near-maximum Taylor expansion 
given in section 2.3. 
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(ix) Use the apparent horizon relations to perform an error correction, as explained in 
section 4. 

(x) Once the data is corrected, start the procedure from the beginning, with the 
adjusted data as input. 

(xi) In practice steps ix-x need to be repeated until the outcome of the algorithm and 
the series expansion match smoothly — see Sec. 4 for details. 

With M (r), E{r), and ^^(r) in hand, the LT model is determined, and one can calculate 
its state and evolution for any instant using (4). 

6. Results 

It should be re-iterated that our present purpose is to vahdate various possible methods 
of reducing observational data to metric information, and investigate how sensitive they 
are to variations in approach. We do not suggest the LT models derived here represent 
reality. Rather we are working towards a reliable method that can be used with much 
larger and more complete survey datasets. The ultimate goal is to have a working 
algorithm that allows the metric of the cosmos to be derived from the observations, 
rather than choosing a metric first and then trying to fit its functions to observations. 

Let us now test different combinations of fitting functions for Tj & Dj (where i 
and j correspond to a fitting type — see Sec. 3). We will consider both models with 
and without the cosmological constant, so all together we have 30 different models. 
The results in form of 0,r{f) Siad tB{r) are presented in Figs. 5 and 6 respectively. 
The results for models with A = are presented using solid lines, results for models 
with the standard model cosmological constant value {Q\ = 0.7) are presented using 
dashed lines. As may be seen, different types of algebraic fits lead to different results 
(scale and range also vary from panel to panel). Clearly, not every pair of r{z) & 
Da{z) fit-functions produces a satisfactory model. Some exhibit shell crossing at the 
current instant (though not along the past null cone): (ri&Z^i), (ri&D2), (T2&-D2), 
(t2 & -D3), and (t3 & D2). Others have M' < for low z which leads to negative density: 
(T2&-D5), (T3&D1), (T3&D2), and (Ts&^Di). While others suffer a similar problem at 
high redshifts: (T1&D5), (t2&L>i), and (T2&r>5). 

6.1. Error estimation 

There are two sources of error in our analysis: the observational errors (both random 
and systematic) , and the (systematic) errors introduced by the choice of function used 
to fit the data. To estimate the 'confidence region', or rather 'error band', we plotted all 
models that do not exhibit a shell crossing, or have M' < 0. In addition we also plotted 
results for models with r ± Ar and Da ± AD^, where At and AD a are uncertainties 
calculated separately for each fit ff. 

tf The error of the quantity Q, that depends on parameters qi, i.e. Q = Q{qi), is estimated using the 
propagation of uncertainty method: (AQ)^ = J2i{'^Q/9Qi)^{^qi)^ ■ 
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Figure 5. Present-day density profile ^^(r) for models with fits to observations in 
the form of Ti{z) and Dj{z) (where i and j correspond to a fitting type — see Sec. 
3). The pair of r and D is given in the upper left corner of each panel. Solid line 
corresponds to a model with Oa = and dashed line to Oa = 0.7. (Note that these 
profiles do not represent observations, which are necessarily on the PNC.) 
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Figure 6. Bang time function tB{r) for models with fits to observations in the form 
of Ti{z) and Dj{z) (where i and j correspond to a fitting type — see Sec. 3). The pair 
of r and D is given in upper or lower left corner of each panel. Solid line corresponds 
to a model with Sl\ = and dashed line to Sl\ = 0.7. 



6.2. Discussion 

The error band and the individual model curves for the present-day density profile are 
presented in Fig. 7, and those for the bang time function in Fig- 8. These show that, 
for models with A = 0, the dispersion in the Qm plot is not as large as in the case of 
Q\ = 0.7. Also, with Q\ = 0, a clear increasing trend is visible up to Rq = 3.5 Gpc at 
least, whereas in the case of Q/^ = 0.7 there is far more uncertainty, and no discernible 
trend. As to the bang time function, both cases allow for decreasing, increasing or a 
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Figure 7. Constraints on the present-day density profile from luminosity distance 
and maximum stellar age data. The curves are the individual fits under different 
assumptions. The green region is the envelope of all fitted models satisfying some 
minimal reasonableness requirements (see text). The region is large because there is 
still a lot of scatter and uncertainty in the data. 




Figure 8. Constraints on the bang time function from luminosity distance and 
maximum stellar age data. The curves are the individual fits, and the green region is 
the envelope of all reasonable models. 



constant ts function up to at least 3.5 Gpc. 

The distance of 3.5 Gpc corresponds to z ^ 1.5, and at this distance the fit to the 
data is no longer reliable; as shown in Fig. 2, for z > 1.5 there is only small number of 
data. Therefore the behaviour of these models at larger distances cannot be considered 
at all accurate, though it is interesting to notice that in all cases (with and without the 
cosmological constant) for i?o > 3.5 Gpc there is a definite decreasing trend. 

7. Conclusions 

We have pursued the broad question of what observations can tell us about the geometry 
of our universe, and to what extent the observations support homogeneity. As the 
cosmological datasets become larger, more accurate, and reach deeper z, we should 
be able to relax our assumptions and derive our spacetime metric in increasing detail, 
verifying to what extent, and on what scale, it is homogeneous on average. While some 
assumptions will always be required, it is important to test different sets of assumptions, 
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and to test theories in different ways. 

Here we have shown how measurements of maximum galactic ages, together with 
SNIa luminosities, over the available range of redshifts, can be used to determine the 
free metric functions of a spherically symmetric cosmology — the Lemaitre-Tolman 
metric. Previous work used galaxy number counts instead of maximum ages. Reasons 
why a spherically symmetric model is a good place to start such an investigation, 
and appropriate to currently available datasets, are discussed above and elsewhere 
[16, 56, 63, 44]. The algorithm does a numerical integration of the DEs we derived, 
but requires careful treatment of the origin, z — 0, and the neighbourhood of the 
maximum in the diameter distance, R = Rm — the apparent horizon, to avoid numerical 
difficulties. We have adapted the apparent horizon correction method used for a different 
dataset combination [43, 56, 63] to our case, and shown how the properties of this 
apparent horizon provide relationships between observables and integrated quantities 
that allow one to detect a cumulative systematic error in the data and make a correction, 
thus ensuring a self-consistent solution. All attempts to integrate past null cone data 
face potential divergences at the apparent horizon. We have discussed the problem of 
extracting smooth functions from discrete data, and estimating the uncertainties in the 
output. We adopted a method, based on a range of possible fitting functions, that is 
appropriate to the limited and noisy data currently available. 

Also, in contrast to the inverse methods usually considered, we did not limit 
ourselves to the zcro-A LT models. Usually, when considering inhomogeneous models, 
people arc seeking to explain the data without the cosmological constant. Recently 
however, inhomogeneous models with the cosmological constant have started to be 
considered more seriously; see for example [61] or [14]. In order to fully test the 
assumption of homogeneity we cannot set A = 0, thus we consider two cases = 
and Qa = 0.7. 

We have executed our algorithm using currently available data, and presented the 
results. Though the data is not yet good enough to attach any confidence to our output, 
we have analysed our graphs as an interesting exercise. In fact our results are quite 
reasonable, and in particular we have shown that there is no conflict between the data 
used and an LT model without A. We expect much stronger constraints to emerge 
as the data improves. It should be emphasised that the stated confidence levels in 
the concordance model are much higher exactly because the homogeneity assumption 
removes huge amounts of potential uncertainty. Homogeneity is not being tested in 
those confidence calculations. We note that comparisons of north and south values in 
e.g. the Sloan data, or angular variations in the Hubble fiow, have yielded differences 
of uncertain significance. 

Thus we have demonstrated that the metric extraction method is both theoretically 
and numerically viable, and should be developed and extended, as an important and 
strong check on homogeneity. It is because both time and distance vary down our past 
null cone in the radial direction, so observations are affected by source evolution as well as 
cosmic expansion, that testing radial homogeneity on top of that — without assuming it 
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somewhere — is very difficult. The relationship between the the paths of light rays — the 
past null cone — and the geometry (metric) they pass through is intrinsically non-linear. 
Therefore our analysis of observations must take this into account. The generalisation 
of the metric extraction method to handling angular variation is intrinsically complex, 
but must be pursued. 

In [54] the authors took an LT model that reproduces the luminosity and density 
curves of the currently preferred ACDM model. They estimated the age of the universe 
at the central location (at us) for this LT model, and with their method they found the 
age of their chosen LT model near us is < 11.7 Gy at a la confidence level, but still 
within la of the locally measured age limits. We point out that they did not find the 
best fit of an LT model to the observations. They assumed the luminosity and density 
curves of the best fit FLRW model are the best fit for any model. However, once one 
allows inhomogeneity, as one should if LT is being used, there is much more leeway 
in the luminosity and density data curves than is apparent from the best-fit FLRW 
parameters and the quoted confidence levels. In contrast, our work uses the data to 
directly construct the model, and finds the age measurements are consistent with a 
range of reasonable inhomogeneous models (including some close to homogeneous). 
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